#include "generalization.h"
#include <iostream>
#include <device_launch_parameters.h>
#include <thrust/device_vector.h>
#include <thrust/sequence.h>
#include <thrust/sort.h>
#include <fstream>
#include <thrust/host_vector.h>
#include <thrust/scan.h>
#include <fstream>
#include <stdlib.h>     
#include <time.h> 
#include <math.h>
#include <set>
#include <sstream>
#include <device_launch_parameters.h>
#include <algorithm>
#include <cuda_profiler_api.h>
using namespace std;
using namespace thrust;


//to record the time spent on those major steps in processing the bi-directional query
vector<double> exec_time;

// compare upper bound of Result object
struct Ubcomapre
{
	__host__ __device__ bool operator()(const Result &lhs, const Result &rhs) const{
		if(lhs.query < rhs.query) return true;
		else if(lhs.query > rhs.query) return false;
		if(lhs.ub == rhs.ub) return lhs.lb < rhs.lb;
		else return lhs.ub < rhs.ub;
	}
};

// compare lower bound of Result object
struct Lbcompare
{
	__host__ __device__ bool operator()(const Result &lhs, const Result &rhs) const{
		if(lhs.query < rhs.query) return true;
		else if(lhs.query > rhs.query) return false;
		return lhs.lb < rhs.lb;
	}
};

void print_seperate_line(){ cout << "++++++++++++++++++++++++++++++++++++++" << endl;}
void print_result_vec(vector<Result>& result_vec){
	cout << "result of output: " << endl;
	for(int i = 0; i < result_vec.size(); i++){
		result_vec[i].print_result_entry();
		cout << endl;
	}
	print_seperate_line();
}

// allocate GPU memory
template <class T>
void init_gpu_memory(host_vector<T>& hvec, device_vector<T>** dvec){
	*dvec = new device_vector<T>(hvec.size());
	**dvec = hvec;
}

// round float to a close integer
__device__ int rounding(float x){return (int)(x+0.5);}

// distance function for values on the same dimension
__device__ T distance_func(float a, float b, int func_type, float bound, float dim_weight){
	float result;
	float diff = fabs(a - b);
	switch(func_type){
	case 1:
	{
		result = diff;
		
	}
	default:
	{
		result = (diff*diff);
	}
	}
	result = (diff < bound) ? 1 : result; // the bound step function
	result = (diff == 0) ? 0 : result; // make sure the 0 case is correct
	result *= dim_weight; // add dimension weight

	return result;
}

__device__ void blk_sort_shared_memory(float* shared_array, float* temp_array, int size){
	int rounds = ceil(log2((double)size));

	int base = 1;
	for(int i = 1; i <= rounds; i++){
		
		if(threadIdx.x % (2*base) == 0){
			int start_idx_x = threadIdx.x;
			int end_idx_x = start_idx_x + base;
			end_idx_x = end_idx_x < size ? end_idx_x : size;

			int start_idx_y = end_idx_x;
			int end_idx_y = start_idx_y + base;
			end_idx_y = end_idx_y < size ? end_idx_y : size;

			int x_pointer = start_idx_x;
			int y_pointer = start_idx_y;

			int output_ptr = x_pointer;
			while(x_pointer < end_idx_x || y_pointer < end_idx_y){
				if(x_pointer >= end_idx_x)
					temp_array[output_ptr++] = shared_array[y_pointer++];
				else if(y_pointer >= end_idx_y)
					temp_array[output_ptr++] = shared_array[x_pointer++];
				else if(shared_array[x_pointer] < shared_array[y_pointer])
					temp_array[output_ptr++] = shared_array[y_pointer++];
				else
					temp_array[output_ptr++] = shared_array[x_pointer++];
			}			
		}
		__syncthreads();
		if(threadIdx.x < size){
			shared_array[threadIdx.x] = temp_array[threadIdx.x];
		}
		__syncthreads();
		base *= 2;
	}
}

// computing the counting while saving the position
__global__ void compute_mapping_saving_pos(QueryInfo* query_list,InvlistEnt* invert_list,int* invert_list_idx,QueryFeatureEnt* query_feature,bool point_search){
	int bid = blockIdx.x; // mapped to different dimension of a query
	int tid = threadIdx.x;
	int block_size = blockDim.x;

	QueryInfo* query = &query_list[0];

	if(bid < query->numOfDimensionToSearch){
	//for(int iter = 0; iter < query->numOfDimensionToSearch; iter++){
		int dim = query->searchDim[bid];
		float dim_weight = query->dimWeight[bid];
		int dim_dist_func = query->distanceFunc[bid];
		float dim_lb_dist_func = query->lowerBoundDist[bid];
		float dim_ub_dist_func = query->upperBoundDist[bid];
		int bound_down_search = query->lowerBoundSearch[bid];
		int bound_up_search = query->upperBoundSearch[bid];

		int down_pos = query->lastPos[dim].x;
		int up_pos = query->lastPos[dim].y;

		int y_dim_value = rounding(query->keyword[dim]);

		int keyword = y_dim_value + dim * MAX_DIM_VALUE;

		if(down_pos == 0 && up_pos == 0){
			int invert_list_start = keyword == 0 ? 0 : invert_list_idx[keyword-1];
			int invert_list_end = invert_list_idx[keyword];
			int invert_list_size = invert_list_end - invert_list_start;
			int process_round = invert_list_size / block_size + (invert_list_size % block_size != 0);

			for(int j = 0; j < process_round; j++){
				int idx = invert_list_start+j*block_size+tid;
				if(idx < invert_list_end){
					InvlistEnt inv_ent = invert_list[idx];
					//int target_idx = bid*invert_list_spec.maxFeatureID+inv_ent;
					//query_feature[target_idx].count += 1;
					atomicAdd(&(query_feature[inv_ent].count),1);
				}
			}
		}
		if(point_search){
			//__syncthreads();
			//continue;
			return;
		}

		// going downward
		int down_compute = invert_list_spec.numOfDocToExpand;
		int down_value = y_dim_value - down_pos; // move the position that last iteartion possessed.

		while(down_value - 1 >= invert_list_spec.minDomainForAllDimension[dim]  // make sure dimension is above minimum dimension value
			&& down_compute > 0	// make sure the number of compute element is above 0
			&& down_value > y_dim_value - bound_down_search
			){
			down_value--;
			int down_keyword = down_value + dim * MAX_DIM_VALUE;

			int invert_list_start = down_keyword == 0 ? 0 : invert_list_idx[down_keyword-1];
			int invert_list_end = invert_list_idx[down_keyword];
			int invert_list_size = invert_list_end - invert_list_start;
			int process_round = invert_list_size / block_size + (invert_list_size % block_size != 0);

			// calcuate the distance from the current point to the query along this dimension
			float true_dist = distance_func(y_dim_value,down_value,dim_dist_func,dim_lb_dist_func,dim_weight);

			for(int j = 0; j < process_round; j++){
				int idx = invert_list_start+j*block_size+tid;
				if(idx < invert_list_end){
					InvlistEnt inv_ent = invert_list[idx];
					//int target_idx = bid*invert_list_spec.maxFeatureID+inv_ent;
					//query_feature[target_idx].count += 1;
					//query_feature[target_idx].ACD += true_dist;
					atomicAdd(&(query_feature[inv_ent].count),1);
					atomicAdd(&(query_feature[inv_ent].ACD),true_dist);
					//if(bid == 1) printf("query %d meet doc %d going down\n",bid,inv_ent);
				}
			}

			down_compute -= invert_list_size;
		}

		// going upward
		int up_compute = invert_list_spec.numOfDocToExpand;
		int up_value = y_dim_value + up_pos; // move the position that last iteartion possessed.

		while(up_value + 1 <= invert_list_spec.maxDomainForAllDimension[dim] // make sure dimension is below maximum dimension value
			&& up_compute > 0 // make sure the number of compute element is above 0
			&& up_value < y_dim_value + bound_up_search
			){
			up_value++;
			int up_keyword = up_value + dim * MAX_DIM_VALUE;

			int invert_list_start = up_keyword == 0 ? 0 : invert_list_idx[up_keyword-1];
			int invert_list_end = invert_list_idx[up_keyword];
			int invert_list_size = invert_list_end - invert_list_start;
			int process_round = invert_list_size / block_size + (invert_list_size % block_size != 0);

			// calcuate the distance from the current point to the query along this dimension
			float true_dist = distance_func(y_dim_value,up_value,dim_dist_func,dim_ub_dist_func,dim_weight);

			for(int j = 0; j < process_round; j++){
				int idx = invert_list_start+j*block_size+tid;
				if(idx < invert_list_end){
					InvlistEnt inv_ent = invert_list[idx];
					/*int target_idx = bid*invert_list_spec.maxFeatureID+inv_ent;
					query_feature[target_idx].count += 1;
					query_feature[target_idx].ACD += true_dist;*/
					atomicAdd(&(query_feature[inv_ent].count),1);
					atomicAdd(&(query_feature[inv_ent].ACD),true_dist);
					//if(bid == 1) printf("query %d meet doc %d going up\n",bid,inv_ent);
				}
			}

			up_compute -= invert_list_size;
		}
		//if(bid == 1 && tid == 0) printf("++++++++++++++++++++++++++++++++++\n");

		__syncthreads();
		if(tid == 0){
			query->lastPos[dim].x =  y_dim_value - down_value;
			query->lastPos[dim].y = up_value - y_dim_value;
		}
	}
}

// count_vec is for the output position for each thread
// blk_count_vec is for the later processing for selecting topK document for each query since each block correspond to a query
__global__ void prefix_count(QueryFeatureEnt* query_feature,int* count_vec, int* blk_count_vec,int threshold){
	int bid = blockIdx.x;
	int tid = threadIdx.x;
	int oid = bid * blockDim.x + tid;

	int block_start = bid * invert_list_spec.maxFeatureID;
	int block_end = block_start + invert_list_spec.maxFeatureID;
	int round = invert_list_spec.maxFeatureID / blockDim.x + (invert_list_spec.maxFeatureID % blockDim.x != 0);


	__shared__ int shared_count;
	if(tid == 0) shared_count = 0;
	__syncthreads();

	int count = 0;
	for(int i = 0; i < round; i++){
		int idx = block_start + i * blockDim.x + tid;
		if(idx < block_end){
			if(query_feature[idx].count > threshold){
				count++;
				atomicAdd(&shared_count,1);
			}
		}
	}
	count_vec[oid] = count;
	__syncthreads();

	if(tid == 0) blk_count_vec[bid] = shared_count;
}

// this function is used for point search query
__global__ void output_result_point_search(
	QueryFeatureEnt* query_feature,
	QueryInfo* query_set, // remeber the counting position
	int* ending_idx,Result* result_vec,
	int threshold){

	int pos = blockIdx.x * blockDim.x + threadIdx.x;
	
	int start_idx = pos == 0 ? 0 : ending_idx[pos-1];

	int block_start = blockIdx.x * invert_list_spec.maxFeatureID;
	int block_end = block_start + invert_list_spec.maxFeatureID;
	int round = invert_list_spec.maxFeatureID / blockDim.x + (invert_list_spec.maxFeatureID % blockDim.x != 0);

	for(int i = 0; i < round; i++){
		int idx = block_start + i * blockDim.x + threadIdx.x;
		if(idx < block_end){
			float count = query_feature[idx].count;
			if(count > threshold){
				Result new_result;
				new_result.query = blockIdx.x;
				new_result.feature_id = i * blockDim.x + threadIdx.x;
				new_result.count = count;
				result_vec[start_idx] = new_result;
				start_idx++;
			}
		}
	}
}

// this function is used for bi directional search query
__global__ void output_result_bidrection_search(
	QueryFeatureEnt* query_feature,
	QueryInfo* query_set, // remeber the counting position
	int* ending_idx,Result* result_vec,
	int threshold){

	int pos = blockIdx.x * blockDim.x + threadIdx.x;
	
	int start_idx = pos == 0 ? 0 : ending_idx[pos-1];

	int block_start = blockIdx.x * invert_list_spec.maxFeatureID;
	int block_end = block_start + invert_list_spec.maxFeatureID;
	int round = invert_list_spec.maxFeatureID / blockDim.x + (invert_list_spec.maxFeatureID % blockDim.x != 0);
	int block_num_search_dim = query_set[blockIdx.x].numOfDimensionToSearch;

	__shared__ float queryLowerBound[MAX_DIM];
	__shared__ float queryUpperBound[MAX_DIM];
	__shared__ float temp_shared[MAX_DIM];

	if(threadIdx.x < block_num_search_dim){
		int dim = query_set[blockIdx.x].searchDim[threadIdx.x];
		int query_dim_value = query_set[blockIdx.x].keyword[dim];
		int dim_dist_func = query_set[blockIdx.x].distanceFunc[threadIdx.x];
		float dim_weight = query_set[blockIdx.x].dimWeight[threadIdx.x];
		float dim_lb_dist_func = query_set[blockIdx.x].lowerBoundDist[threadIdx.x];
		float dim_ub_dist_func = query_set[blockIdx.x].upperBoundDist[threadIdx.x];
		float min_bound_dist_func;
		float max_bound_dist_func;

		// make sure the min distance function bound is lower and max distance function is higher
		if(dim_lb_dist_func > dim_ub_dist_func){
			min_bound_dist_func = dim_ub_dist_func;
			max_bound_dist_func = dim_lb_dist_func;
		}else{
			min_bound_dist_func = dim_lb_dist_func;
			max_bound_dist_func = dim_ub_dist_func;
		}

		int2 q_pos = query_set[blockIdx.x].lastPos[dim];
		// make sure the bound is correct when upward and downward search all reach the maximum
		bool all_reach_max = (q_pos.x == query_dim_value - invert_list_spec.minDomainForAllDimension[dim]) && 
			(q_pos.y == (invert_list_spec.maxDomainForAllDimension[dim] - query_dim_value));


		int min_lb = all_reach_max ? max(q_pos.x,q_pos.y) : min(q_pos.x,q_pos.y); // modified to the min compare to modified.cu which use max
		int min_value = 0;
		queryLowerBound[threadIdx.x] = distance_func(min_lb,min_value,dim_dist_func,min_bound_dist_func,dim_weight);

		
		int max_ub = max(query_dim_value - invert_list_spec.minDomainForAllDimension[dim],invert_list_spec.maxDomainForAllDimension[dim]-query_dim_value);
		queryUpperBound[threadIdx.x] = distance_func(max_ub,min_value,dim_dist_func,max_bound_dist_func,dim_weight);
	}

	__syncthreads();

	// sort the lower and up bound
	blk_sort_shared_memory(queryLowerBound,temp_shared,invert_list_spec.totalDimension);
	blk_sort_shared_memory(queryUpperBound,temp_shared,invert_list_spec.totalDimension);
	
	// calculate lower bound for this query, total distance for current retrieve point
	int lower_bound_sum = 0;
	for(int i = 0; i < block_num_search_dim; i++) lower_bound_sum += queryLowerBound[i];

	for(int i = 0; i < round; i++){
		int idx = block_start + i * blockDim.x + threadIdx.x;
		if(idx < block_end){
			int count = query_feature[idx].count;
			if(count > threshold){
				Result new_result;
				new_result.query = blockIdx.x;
				new_result.feature_id = i * blockDim.x + threadIdx.x;
				new_result.count = count;
				
				float ACD = query_feature[idx].ACD;

				// DL(f) = DL(us) - sum(max count number of dt) + ACD(f)
				float lb = lower_bound_sum;
				for(int m = 0; m < count; m++){
					lb -= queryLowerBound[m];
				}
				lb += ACD;

				// DU(f) = ACD(f) + sum(max 128 - c of upper bound)
				float ub = ACD;
				for(int m = 0; m < block_num_search_dim - count; m++){
					ub += queryUpperBound[m];
				}

				new_result.lb = lb;
				new_result.ub = ub;

				result_vec[start_idx] = new_result;
				start_idx++;
			}
		}
	}
}

__global__ void terminate_check(Result* ub_sorted,Result* lb_sorted, int* end_idx, int* output, int K, int round_num){
	int tid = threadIdx.x;
	int bid = blockIdx.x;

	int blk_start_idx = bid == 0 ? 0 : end_idx[bid-1];
	int blk_end_idx = end_idx[bid];

	if(blk_end_idx - blk_start_idx <= K) return;

	blk_end_idx = blk_start_idx + K + 1;

	int round = (K+1) / blockDim.x + ((K+1) % blockDim.x != 0);

	__shared__ int termination_check;
	if(tid == 0) termination_check = 0;
	__syncthreads();

	for(int i = 0; i < round; i++){
		int idx = blk_start_idx + i * blockDim.x + tid;
		if(idx < blk_end_idx){
			Result lb_elem = lb_sorted[idx];
			bool duplicated = false; // if element in the lb sorted array existed in the ub sorted array
			for(int j = 0; j < K; j++){
				Result ub_elem = ub_sorted[blk_start_idx + j];
				if(lb_elem.feature_id == ub_elem.feature_id){
					duplicated = true;
					break;
				}
			}
			Result ub_elem = ub_sorted[blk_start_idx + K - 1];
			bool lb_greater_ub = (lb_elem.lb >= ub_elem.ub);

			if(!duplicated && lb_greater_ub){
				//printf("reached here tid %d and bid %d\n",tid,bid);
				atomicAdd(&termination_check,1);
			}
		}
	}
	__syncthreads();
	if(tid == 0){
		if(termination_check != 0){
			output[bid] = 0;
		}
	}
}

__global__ void extract_topK(Result* ub_sorted,int* end_idx,Result* output, int K){
	int tid = threadIdx.x;
	int bid = blockIdx.x;
	int oid = bid * K + tid;

	int blk_start_idx = bid == 0 ? 0 : end_idx[bid-1];
	int blk_end_idx = blk_start_idx + K;

	int round = K / blockDim.x + (K % blockDim.x != 0);

	for(int i = 0; i < round; i++){
		int idx = blk_start_idx + i * blockDim.x + tid;
		if(idx < blk_end_idx){
			output[oid] = ub_sorted[idx];
			oid += blockDim.x;
		}
	}
}

QueryManager::QueryManager(QuerySpecification& query_spec,int rand_invert_list_size){
	srand(1);
	invert_list_spec_host.totalDimension = query_spec.totalDimension;
	if(invert_list_spec_host.totalDimension <= 0 || invert_list_spec_host.totalDimension > MAX_DIM){
		cout << "total number of dimension of this invert list is wrong" << endl;
		exit(1);
	}
	invert_list_spec_host.numOfDocToExpand = query_spec.numOfDocToExpand;
	if(invert_list_spec_host.numOfDocToExpand <= 0){
		cout << "cannot handle numOfDocToExpand is equal or below 0" << endl;
		exit(1);
	}
	for(int i = 0; i < invert_list_spec_host.totalDimension; i++){
		invert_list_spec_host.minDomainForAllDimension[i] = query_spec.minDomainForAllDimension[i];
		invert_list_spec_host.maxDomainForAllDimension[i] = query_spec.maxDomainForAllDimension[i];
		if(invert_list_spec_host.minDomainForAllDimension[i] < 0 || invert_list_spec_host.maxDomainForAllDimension[i] >= MAX_DIM_VALUE ||
			invert_list_spec_host.minDomainForAllDimension[i] >= invert_list_spec_host.maxDomainForAllDimension[i]){
			cout << "max and min value of dimension " << i << " has problem: should be positive and under MAX_DIM_VALUE" << endl;
			exit(1);
		}
	}

	if(query_spec.invertedListPath.empty() && query_spec.invertedIndexLengthFile.empty()){
		invert_list_spec_host.maxFeatureID = rand_invert_list_size;
		rand_inv_list();
	}else{
		// read inverted list from file
		//this->get_invert_list_from_file(query_spec.invertedListPath,query_spec.invertedIndexLengthFile);
		get_invert_list_from_binary_file(query_spec.invertedListPath);
	}

	this->d_query_feature = NULL;
	this->d_query_info = NULL;
	this->d_valid_query = NULL;

	this->topK = query_spec.topK;
	if(topK <= 0){
		cout << "topK should be a positive number" << endl;
		exit(1);
	}
	this->query_initialized = false;
}

void QueryManager::rand_inv_list(){
	vector<vector<int> > keyword_feature(invert_list_spec_host.totalDimension * MAX_DIM_VALUE);
	for(int i = 0; i < invert_list_spec_host.maxFeatureID; i++){
		for(int x = 0; x < invert_list_spec_host.totalDimension;  x++){
			int dim_range = invert_list_spec_host.maxDomainForAllDimension[x] - invert_list_spec_host.minDomainForAllDimension[x];
			int dim_base = invert_list_spec_host.minDomainForAllDimension[x];
			int dim_value = rand() % dim_range + dim_base;
			int keyword = x * MAX_DIM_VALUE + dim_value;
			keyword_feature[keyword].push_back(i);
		}
	}

	host_vector<InvlistEnt> h_invert_list;
	host_vector<int> h_invert_list_idx;
	int count = 0;
	// construct inverted list
	for(int i = 0; i < invert_list_spec_host.totalDimension * MAX_DIM_VALUE; i++){
		for(int j = 0; j < keyword_feature[i].size(); j++){
			InvlistEnt new_ent;
			new_ent = keyword_feature[i][j];
			h_invert_list.push_back(new_ent);
		}
		count += keyword_feature[i].size();
		h_invert_list_idx.push_back(count);
	}
	
	// transfer to gpu
	init_gpu_memory(h_invert_list,&d_invert_list);
	init_gpu_memory(h_invert_list_idx,&d_invert_list_idx);
}

void QueryManager::rand_query(int query_num){
	host_vector<QueryInfo> h_query_info(query_num);
	for(int i = 0; i < query_num; i++){
		// setup search dimension
		/*h_query_info[i].numOfDimensionToSearch = rand() % invert_list_spec_host.totalDimension + 1;
		set<int> temp_set;
		while(temp_set.size() != h_query_info[i].numOfDimensionToSearch)
			temp_set.insert(rand() % invert_list_spec_host.totalDimension);
		int temp_pos = 0;
		for(set<int>::iterator it = temp_set.begin(); it != temp_set.end(); it++)
			h_query_info[i].searchDim[temp_pos++] = *it;*/
		h_query_info[i].numOfDimensionToSearch = invert_list_spec_host.totalDimension;
		for(int j = 0; j < invert_list_spec_host.totalDimension; j++) 
			h_query_info[i].searchDim[j] = j;

		// setup query keyword
		for(int x = 0; x < invert_list_spec_host.totalDimension; x++){
			int dim_range = invert_list_spec_host.maxDomainForAllDimension[x] - invert_list_spec_host.minDomainForAllDimension[x];
			int dim_base = invert_list_spec_host.minDomainForAllDimension[x];
			float dim_value = rand() % dim_range + dim_base;
			h_query_info[i].keyword[x] = dim_value;
		}
	}
	init_gpu_memory(h_query_info,&d_query_info);
}

void QueryManager::init_new_query(vector<GpuQuery>& query_set,int rand_query_size){
	clear_query_memory();
	
	if(query_set.empty()){
		// TO BE DEVELOPED: has not setup query dimWeight, distanceFunc, upperBoundDist and lowerBoundDist
		// they are using default values now
		this->invert_list_spec_host.numOfQuery = rand_query_size;
		rand_query(this->invert_list_spec_host.numOfQuery);
	}else{
		this->invert_list_spec_host.numOfQuery = query_set.size();
		// read from query structure
		host_vector<QueryInfo> h_query_info(this->invert_list_spec_host.numOfQuery);
		for(int i = 0; i < this->invert_list_spec_host.numOfQuery; i++){
			// copy keywords
			for(int j = 0; j < this->invert_list_spec_host.totalDimension; j++)
				h_query_info[i].keyword[j] = query_set[i].keywords[j];

			// set up dimension
			h_query_info[i].numOfDimensionToSearch = query_set[i].numOfDimensionToSearch;
			h_query_info[i].aggregateFunc = query_set[i].aggregateFunction;
			for(int j = 0; j < h_query_info[i].numOfDimensionToSearch; j++){
				h_query_info[i].searchDim[j] = query_set[i].dimensionSet[j].dimension;
				h_query_info[i].dimWeight[j] = query_set[i].dimensionSet[j].weight;
				h_query_info[i].distanceFunc[j] = query_set[i].dimensionSet[j].distanceFunctionPerDimension;
				
				// MAY NEED TO CHANGE LATER: since each search dim may have different bound
				h_query_info[i].upperBoundDist[j] = query_set[i].upwardSearchBound;
				h_query_info[i].lowerBoundDist[j] = query_set[i].downwardSearchBound;

				// NEED TO ADD LATER: search lower bound and upper bound of different dimension.
			}
		}
		init_gpu_memory(h_query_info,&d_query_info);
	}

	cudaMemcpyToSymbol(invert_list_spec,&this->invert_list_spec_host,sizeof(InvertListSpecGPU));

	// checking memory size
	if((double)invert_list_spec_host.maxFeatureID*invert_list_spec_host.numOfQuery*sizeof(QueryFeatureEnt) > std::numeric_limits<int>::max()){
		cout << "too much memory is used for query feature mapping" << endl;
		exit(1);
	}

	d_query_feature = new device_vector<QueryFeatureEnt>((unsigned int)invert_list_spec_host.maxFeatureID*invert_list_spec_host.numOfQuery);
	d_valid_query = new device_vector<int>(invert_list_spec_host.numOfQuery,1);
	num_of_rounds = 0;

	this->query_initialized = true;
	cout << "finish initialization" << endl;
}

QueryManager::~QueryManager(){
	delete d_invert_list;
	delete d_invert_list_idx;
	clear_query_memory();
}

void QueryManager::clear_query_memory(){
	if(d_query_info != NULL) delete d_query_info;
	if(d_query_feature != NULL) delete d_query_feature;
	if(d_valid_query != NULL) delete d_valid_query;
}

void QueryManager::print_query(){
	cout << "query info" << endl;
	host_vector<QueryInfo> host_query_info = *(this->d_query_info);
	for(int i = 0; i < host_query_info.size(); i++){
		cout << "query " << i << ":";
		for(int j = 0; j < this->invert_list_spec_host.totalDimension; j++)
			cout << host_query_info[i].keyword[j] << " ";
		cout << endl;
	}
	print_seperate_line();
}

void QueryManager::print_query_doucment_mapping(){
	cout << "query document mapping" << endl;
	host_vector<QueryFeatureEnt> host_mapping = *(this->d_query_feature);
	for(int i = 0; i < invert_list_spec_host.numOfQuery; i++){
		for(int j = 0; j < invert_list_spec_host.maxFeatureID; j++){
			QueryFeatureEnt entry = host_mapping[i*invert_list_spec_host.maxFeatureID+j];
			cout << entry.count << " ";
		}
		cout << endl;
	}
	print_seperate_line();
}

void QueryManager::print_invert_list(){
	cout << "invert list" << endl;
	host_vector<InvlistEnt> host_invert_list = *(this->d_invert_list);
	host_vector<InvlistEnt> host_invert_list_idx = *(this->d_invert_list_idx);
	for(int i = 0; i < host_invert_list_idx.size(); i++){
		int dim = i / MAX_DIM_VALUE;
		int value = i % MAX_DIM_VALUE;
		cout << "(" << dim << "," << value << "):";
		int start = i == 0 ? 0 : host_invert_list_idx[i-1];
		int end = host_invert_list_idx[i];
		for(int j = start; j < end; j++) cout << host_invert_list[j] << " ";
		cout << endl;
	}
	print_seperate_line();
}

bool QueryManager::bi_direction_query(vector<Result>& result,int threshold){
	if(!this->query_initialized){
		cout << "query has not been initialized" << endl;
		exit(1);
	}
	if(THREAD_PER_BLK < MAX_DIM){
		cout << "thread per block must be larger than max dimension of x" << endl;
		exit(1);
	}
	if(this->invert_list_spec_host.numOfQuery > 65535){
		cout << "number of query should not exceed 65535: need to modify the code" << endl;
		exit(1);
	}
	if(THREAD_PER_BLK > 1024){
		cout << "THREAD_PER_BLK should not exceed 1024: need to modify the code" << endl;
		exit(1);
	}

	
	cudaProfilerStart();			
	cudaEvent_t start, stop;
	float elapsedTime;
	cudaEventCreate(&start);
	cudaEventCreate(&stop);

	cudaEventRecord(start, 0); // ******************** //  timing section
	
	// assume there is only one query
	compute_mapping_saving_pos<<<MAX_DIM,THREAD_PER_BLK>>>(
		raw_pointer_cast(this->d_query_info->data()),
		raw_pointer_cast(this->d_invert_list->data()),
		raw_pointer_cast(this->d_invert_list_idx->data()),
		raw_pointer_cast(this->d_query_feature->data()),false);
	
	device_vector<int> threshold_count(invert_list_spec_host.numOfQuery*THREAD_PER_BLK,0);	// prefix count for each thread
	device_vector<int> query_result_count(invert_list_spec_host.numOfQuery,0);				// prefix count for each block (i.e. each query)

	cudaEventRecord(stop, 0);
	cudaEventSynchronize(stop);
	cudaEventElapsedTime(&elapsedTime, start, stop); 
	exec_time[0] += elapsedTime;

	cudaEventRecord(start, 0); // ******************** //  timing section

	prefix_count<<<this->invert_list_spec_host.numOfQuery,THREAD_PER_BLK>>>(
		raw_pointer_cast(this->d_query_feature->data()),
		raw_pointer_cast(threshold_count.data()),
		raw_pointer_cast(query_result_count.data()),threshold);

	thrust::inclusive_scan(threshold_count.begin(),threshold_count.end(),threshold_count.begin()); // per thread inclusive scan
	thrust::inclusive_scan(query_result_count.begin(),query_result_count.end(),query_result_count.begin()); // per block inclusive scan

	cudaEventRecord(stop, 0);
	cudaEventSynchronize(stop);
	cudaEventElapsedTime(&elapsedTime, start, stop); 
	exec_time[1] += elapsedTime;   // ******************** //
	
	
		
	device_vector<Result> d_result_ub_sorted(threshold_count[threshold_count.size()-1]);

	cudaEventRecord(start, 0);  // ******************** //
	
	output_result_bidrection_search<<<this->invert_list_spec_host.numOfQuery,THREAD_PER_BLK>>>(
		raw_pointer_cast(this->d_query_feature->data()),
		raw_pointer_cast(this->d_query_info->data()),
		raw_pointer_cast(threshold_count.data()),
		raw_pointer_cast(d_result_ub_sorted.data()),
		threshold);

	cudaEventRecord(stop, 0);
	cudaEventSynchronize(stop);
	cudaEventElapsedTime(&elapsedTime, start, stop); 
	exec_time[2] += elapsedTime;
	
	device_vector<Result> d_result_lb_sorted(d_result_ub_sorted.begin(),d_result_ub_sorted.end());

	cudaEventRecord(start, 0);  // ******************** //
	
	thrust::sort(d_result_ub_sorted.begin(),d_result_ub_sorted.end(),Ubcomapre());
	thrust::sort(d_result_lb_sorted.begin(),d_result_lb_sorted.end(),Lbcompare());

	cudaEventRecord(stop, 0);
	cudaEventSynchronize(stop);
	cudaEventElapsedTime(&elapsedTime, start, stop); 
	exec_time[3] += elapsedTime;
	
	
	cudaEventRecord(start, 0);  // ******************** //
	
	terminate_check<<<this->invert_list_spec_host.numOfQuery,THREAD_PER_BLK>>>(
		raw_pointer_cast(d_result_ub_sorted.data()),
		raw_pointer_cast(d_result_lb_sorted.data()),
		raw_pointer_cast(query_result_count.data()),
		raw_pointer_cast(d_valid_query->data()),this->topK,this->num_of_rounds);

	int terminate_sum = thrust::reduce(d_valid_query->begin(),d_valid_query->end());

	cudaEventRecord(stop, 0);
	cudaEventSynchronize(stop);
	cudaEventElapsedTime(&elapsedTime, start, stop); 
	exec_time[4] += elapsedTime;
	cudaEventDestroy(start);
	cudaEventDestroy(stop);
	
	cudaProfilerStop();
	
	device_vector<Result> d_result(this->invert_list_spec_host.numOfQuery*this->topK);

	extract_topK<<<this->invert_list_spec_host.numOfQuery,THREAD_PER_BLK>>>(
		raw_pointer_cast(d_result_ub_sorted.data()),
		raw_pointer_cast(query_result_count.data()),
		raw_pointer_cast(d_result.data()),this->topK);

	num_of_rounds++;

	result.resize(d_result.size());
	cudaMemcpy(&result[0],raw_pointer_cast(d_result.data()),d_result.size()*sizeof(Result),cudaMemcpyDeviceToHost);

	// no more query should be made to this batch
	if(terminate_sum == 0) this->query_initialized = false;

	return terminate_sum == 0;
}

void QueryManager::point_query(vector<Result>& result,int threshold){
	if(!this->query_initialized){
		cout << "query has not been initialized" << endl;
		exit(1);
	}
	if(THREAD_PER_BLK < MAX_DIM){
		cout << "thread per block must be larger than max dimension of x" << endl;
		exit(1);
	}
	if(this->invert_list_spec_host.numOfQuery > 65535){
		cout << "number of query should not exceed 65535: need to modify the code" << endl;
		exit(1);
	}
	if(THREAD_PER_BLK > 1024){
		cout << "THREAD_PER_BLK should not exceed 1024: need to modify the code" << endl;
		exit(1);
	}
	
	compute_mapping_saving_pos<<<this->invert_list_spec_host.numOfQuery,THREAD_PER_BLK>>>(
		raw_pointer_cast(this->d_query_info->data()),
		raw_pointer_cast(this->d_invert_list->data()),
		raw_pointer_cast(this->d_invert_list_idx->data()),
		raw_pointer_cast(this->d_query_feature->data()),true);

	device_vector<int> threshold_count(invert_list_spec_host.numOfQuery*THREAD_PER_BLK,0);	// prefix count for each thread
	device_vector<int> query_result_count(invert_list_spec_host.numOfQuery,0);				// prefix count for each block (query)

	prefix_count<<<this->invert_list_spec_host.numOfQuery,THREAD_PER_BLK>>>(
		raw_pointer_cast(this->d_query_feature->data()),
		raw_pointer_cast(threshold_count.data()),
		raw_pointer_cast(query_result_count.data()),threshold);

	thrust::inclusive_scan(threshold_count.begin(),threshold_count.end(),threshold_count.begin()); // per thread inclusive scan
	thrust::inclusive_scan(query_result_count.begin(),query_result_count.end(),query_result_count.begin()); // per block inclusive scan

	device_vector<Result> d_result(threshold_count[threshold_count.size()-1]);

	output_result_point_search<<<this->invert_list_spec_host.numOfQuery,THREAD_PER_BLK>>>(
		raw_pointer_cast(this->d_query_feature->data()),
		raw_pointer_cast(this->d_query_info->data()),
		raw_pointer_cast(threshold_count.data()),
		raw_pointer_cast(d_result.data()),
		threshold);

	result.resize(d_result.size());
	cudaMemcpy(&result[0],raw_pointer_cast(d_result.data()),d_result.size()*sizeof(Result),cudaMemcpyDeviceToHost);

	this->query_initialized = false; // make sure no other query is called
}

void QueryManager::get_invert_list_from_file(string invert_list_dir,string invert_list_idx_file)
{

	string sorted_file_name(invert_list_idx_file);
	ifstream sorted_file(sorted_file_name.c_str());
	if(!sorted_file.is_open()){
		cout << "sorted index file not found" << endl;
		exit(1);
	}

	cout << "start to read invert list" << endl;

	host_vector<InvlistEnt> h_inv_list;
	host_vector<int> h_inv_list_idx;
	int count = 0;
	int keyword_idx = 0;

	this->invert_list_spec_host.maxFeatureID = -1;

	while(!sorted_file.eof()){
		string input_line;
		getline(sorted_file,input_line);
		stringstream ss(input_line);

		string name;
		int file_idx;
		int file_len;

		ss >> name;
		if(name.empty()) continue;
		ss >> file_len;

		file_idx = atoi(name.c_str());
		string file_name(invert_list_dir);
		file_name.append("/");
		file_name.append(name.c_str());

		int dim_value = file_idx % MAX_DIM_VALUE;
		int dim = file_idx / MAX_DIM_VALUE;
		if(dim_value < invert_list_spec_host.minDomainForAllDimension[dim] ||
			dim_value > invert_list_spec_host.maxDomainForAllDimension[dim]){
			cout << "input invert list has problem for dimension values, should be between min and max dimension value" << endl;
			cout << "dimension " << dim << " min: " << invert_list_spec_host.minDomainForAllDimension[dim] << " and max: " << invert_list_spec_host.maxDomainForAllDimension[dim] << endl;
			cout << "value is: " << dim_value << endl;
			exit(1);
		}

		if(file_idx > invert_list_spec_host.totalDimension * MAX_DIM_VALUE){
			cout << "input invert list does not fit the current configuration, look for MAX_DIM and MAX_DIM_VALUE" << endl;
			sorted_file.close();
			exit(1);
		}

		for(int i = keyword_idx; i < file_idx; i++){
			h_inv_list_idx.push_back(count);
		}

		keyword_idx = file_idx;
		count += file_len;

		ifstream keyword_file(file_name.c_str());
		if(!keyword_file.is_open()){
			cout << "invert list file not found" << endl;
			exit(1);
		}

		for(int i = 0; i < file_len; i++){
			InvlistEnt new_ent;
			keyword_file >> new_ent;
			h_inv_list.push_back(new_ent);
			if(new_ent > this->invert_list_spec_host.maxFeatureID ){
				this->invert_list_spec_host.maxFeatureID  = new_ent;
			}
		}

		keyword_file.close();

	} 
	sorted_file.close();
	cout << "finish reading invert list" << endl;

	this->invert_list_spec_host.maxFeatureID++;

	for(int i = keyword_idx; i < invert_list_spec_host.totalDimension * MAX_DIM_VALUE; i++){
		h_inv_list_idx.push_back(count);
	}

	init_gpu_memory(h_inv_list,&d_invert_list);
	init_gpu_memory(h_inv_list_idx,&d_invert_list_idx);
}


/*
void QueryManager::get_invert_list_from_binary_file(string filename)
{
	
	// initialize the previousKey to be 0;
		//int previousKey = 0;
		int previousKey = 0;
		int key;
		int numberOfFeatures;
		int keyPosition = 0;						
				
		vector<InvlistEnt>  featureIdList;
		vector<int>  featureIdListIndex;				
		
		clock_t start = clock();							  
		ifstream inFile(filename.c_str(), ios::in | ios::binary );
		if ( ! inFile )
		{
			cerr << "error : unable to open input file: " << inFile << endl;
		}
		else if ( ! inFile.is_open() )
		{
			cerr << "error : unable to open input file: " << inFile << endl;
		}
		else
		{
			
			int maxValuePerDim;
			inFile.read((char*) (&maxValuePerDim), sizeof(int) );           
			 
			int numberOfDim;
			inFile.read((char*) (&numberOfDim), sizeof(int) );                                                         
			             
			cout << "numOfDim : " << numberOfDim << endl;
			cout << "maxValuePerDim : " << maxValuePerDim << endl;
			
			// read the key, if reach the end of the file, jump out
			while ( inFile.read((char*) (&key), sizeof(int) ) )
			{			
				int dim_value = key % MAX_DIM_VALUE;
				int dim = key / MAX_DIM_VALUE;
				if(dim > invert_list_spec_host.totalDimension){
						cout << "dim: " << dim << endl;
						cout << "key: " << key << endl;
						cout << "read file problem" << endl;
				}
				if(dim_value < invert_list_spec_host.minDomainForAllDimension[dim] ||
					dim_value > invert_list_spec_host.maxDomainForAllDimension[dim]){
						cout << "input invert list has problem for dimension values, should be between min and max dimension value" << endl;
						cout << "dimension " << dim << " min: " << invert_list_spec_host.minDomainForAllDimension[dim] << " and max: " << invert_list_spec_host.maxDomainForAllDimension[dim] << endl;
						cout << "value is: " << dim_value << endl;
						cout << "key: " << key << endl;
						exit(1);
				}
				
				// add the position for empty keys and the current key
				for ( int j = previousKey; j < key; j++ )
				{
					featureIdListIndex.push_back(keyPosition);
				}

				// read the number of features
				inFile.read((char*) (&numberOfFeatures), sizeof(int) );
				
				// update the current position, and the previousKey
				keyPosition += numberOfFeatures;
				previousKey = key;

				//cout<< "key : " << key <<endl;
				// read the feature ids, and fill in the featureIdList
				if ( numberOfFeatures > 0 )
				{
					int featureID;			
					for ( int j = 0; j < numberOfFeatures; j++ )
					{
						// read the feature ID list
						inFile.read((char*) (&featureID), sizeof(int) );																	
						InvlistEnt new_ent;
						new_ent.feature_id = featureID;
						if(new_ent.feature_id > invert_list_spec_host.maxFeatureID){
							invert_list_spec_host.maxFeatureID = new_ent.feature_id;
						}
						featureIdList.push_back(new_ent);
					}
				}								
			}
			
			MAX_FEATURE_ID++;  //very important!
			
			// fill the key positions for the empty keys
			
			for ( int j = previousKey ; j < MAX_DIM_X * MAX_DIM_Y; j++)
			{
				featureIdListIndex.push_back(keyPosition);
			}

			inFile.close();
			inFile.clear();
			
		}
		
		double init_time = (double)(clock() - start) / CLOCKS_PER_SEC;
		cout<< "Reading Inverted Lists from disk to Host Memory takes: " << init_time << " seconds. " <<endl;
				
		start = clock();
		host_vector<InvlistEnt> featureIdList_h(featureIdList); 
		host_vector<int>  featureIdListIndex_h(featureIdListIndex);
		cout << "size of list: " << featureIdList_h.size() * sizeof(InvlistEnt) << endl;
		cout<< "host vector allocated" <<endl;
		init_gpu_memory(featureIdList_h,&d_invert_list);
		init_gpu_memory(featureIdListIndex_h,&d_invert_list_idx);
		init_time = (double)(clock() - start) / CLOCKS_PER_SEC;
		cout<< "Loading from CPU to GPU takes: " << init_time << " seconds. " <<endl;	
}
*/

//Upated on Jan 13!
void QueryManager::get_invert_list_from_binary_file(string filename){
		// initialize the previousKey to be 0;
		int previousKey = 0;
		int key;
		int numberOfFeatures;
		int keyPosition = 0;
		int offset = 0;
		this->invert_list_spec_host.maxFeatureID = -1;
		int maxNumberOfList = this->invert_list_spec_host.totalDimension * MAX_DIM_VALUE;			
		
		vector<InvlistEnt>  featureIdList;
		vector<int>  featureIdListIndex;
		
		clock_t start = clock();							  
		ifstream inFile(filename.c_str(), ios::in | ios::binary );

		cout << "file name is: " << filename << endl;
		if ( ! inFile )
		{
			cerr << "error : unable to open input file 1: " << inFile << endl;
		}
		else if ( ! inFile.is_open() )
		{
			cerr << "error : unable to open input file 2: " << inFile << endl;
		}
		else
		{
			/*
			//add on Jan 13
			int maxValuePerDim;
			inFile.read((char*) (&maxValuePerDim), sizeof(int) );           
						 
			int numberOfDim;
			inFile.read((char*) (&numberOfDim), sizeof(int) );                                                         
						             
			cout << "numOfDim : " << numberOfDim << endl;
			cout << "maxValuePerDim : " << maxValuePerDim << endl;
			*/
			// read the key, if reach the end of the file, jump out
			while ( inFile.read((char*) (&key), sizeof(int) ) )
			{

				int dim_value = key % MAX_DIM_VALUE;
				int dim = key / MAX_DIM_VALUE;
				if(dim > invert_list_spec_host.totalDimension){
					cout << "dim: " << dim << endl;
					cout << "key: " << key << endl;
					cout << "read file problem" << endl;
				}
				if(dim_value < invert_list_spec_host.minDomainForAllDimension[dim] ||
					dim_value > invert_list_spec_host.maxDomainForAllDimension[dim]){
						cout << "input invert list has problem for dimension values, should be between min and max dimension value" << endl;
						cout << "dimension " << dim << " min: " << invert_list_spec_host.minDomainForAllDimension[dim] << " and max: " << invert_list_spec_host.maxDomainForAllDimension[dim] << endl;
						cout << "value is: " << dim_value << endl;
						cout << "key: " << key << endl;
						exit(1);
				}

				// read the number of features
				inFile.read((char*) (&numberOfFeatures), sizeof(int) );

				// read the feature ids, and fill in the featureIdList
				if ( numberOfFeatures > 0 )
				{
					int featureID;
					for ( int j = 0; j < numberOfFeatures; j++ )
					{
						// read the feature ID list
						inFile.read((char*) (&featureID), sizeof(int) );
						InvlistEnt new_ent = featureID;
						if(new_ent > invert_list_spec_host.maxFeatureID ){
							invert_list_spec_host.maxFeatureID  = new_ent;
						}
						featureIdList.push_back(new_ent);
					}
				}

				// get the offset between keys
				offset = key - previousKey;

				// add the position for empty keys and the current key
				for ( int j = 0; j < offset; j++ )
				{
					featureIdListIndex.push_back(keyPosition);
				}

				// update the current position, and the previousKey
				keyPosition += numberOfFeatures;
				previousKey = key;
			}

			// fill the key positions for the empty keys
			offset = maxNumberOfList - key - 1;
			for ( int j = 0; j < offset; j++ )
			{
				featureIdListIndex.push_back(keyPosition);
			}

			inFile.close();
			inFile.clear();
		}
		
		double init_time = (double)(clock() - start) / CLOCKS_PER_SEC;
		cout<< "Reading Inverted Lists from disk to Host Memory takes: " << init_time << " seconds. " <<endl;
		invert_list_spec_host.maxFeatureID++;

		start = clock();
		host_vector<InvlistEnt> featureIdList_h(featureIdList); 
		host_vector<int>  featureIdListIndex_h(featureIdListIndex);

		cout << "size of list: " << featureIdList_h.size() * sizeof(InvlistEnt) << endl;
		cout<< "host vector allocated" <<endl;
		init_gpu_memory(featureIdList_h,&d_invert_list);
		init_gpu_memory(featureIdListIndex_h,&d_invert_list_idx);
		init_time = (double)(clock() - start) / CLOCKS_PER_SEC;
		cout<< "Loading from CPU to GPU takes: " << init_time << " seconds. " <<endl;
	
}

int main(){
	QuerySpecification test_spec;
	test_spec.numOfDocToExpand = 1000;
	test_spec.topK = 4;
	test_spec.totalDimension = 128;
	
	//test_spec.invertedListPath = "../data/sift_exp_data/sift_small";
	//test_spec.invertedIndexLengthFile = "../data/sift_exp_data/sorted.txt";
	
	test_spec.invertedListPath = "../data/sift_4m/inverted_index.bin";	
	test_spec.maxDomainForAllDimension.resize(test_spec.totalDimension);
	test_spec.minDomainForAllDimension.resize(test_spec.totalDimension);
	for(int i = 0; i < test_spec.totalDimension; i++){
		test_spec.minDomainForAllDimension[i] = 0;            // dangerous:!!
		test_spec.maxDomainForAllDimension[i] = 255;		  // dangerous:!!
	}
	
	clock_t start = clock();
	
	vector<GpuQuery> test_query;

	QueryManager test_manager(test_spec,0);
	test_manager.init_new_query(test_query,1);	

	//test_manager.print_invert_list();
	//test_manager.print_query();

	
	exec_time.resize(5,0);
	double total_time = 0;
	
	vector<Result> temp_result;
	bool finished = false;
	while(!finished){
		finished = test_manager.bi_direction_query(temp_result,0);
		
		double time = (double)(clock() - start) / CLOCKS_PER_SEC;
		total_time += time;
		cout << "iteration takes time: " << time << endl;
		
		print_result_vec(temp_result);
		print_seperate_line();
	}
	
	cout << "finished with total time : " << total_time << endl;
	
	
	//part 0: expansion, update count and ACD; part 1: compute lb and ub for features whose count > 0; part 2: k candidates selection (by sorting in current implementation); part 3: terminate condition check 
	cout << "Profiling the time : " << endl;
		for(int i = 0; i < exec_time.size(); i++){
			cout << "part " << i << " takes " << exec_time[i] << endl;
		}
	
	//vector<Result> temp_result;
	//test_manager.point_query(temp_result,5);
	//print_result_vec(temp_result);
	
	//test_manager.print_query_doucment_mapping();
	
}